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Abstract 

Solvation free energy is an important quantity in Computational Chemistry with a 
variety of applications, especially in drug discovery and design. The accurate predic- 
tion of solvation free energies of small molecules in water is still a largely unsolved 
problem, which is mainly due to the complex nature of the water-solute interac- 
tions. In this letter we develop a scheme for the determination of the electrostatic 
contribution to the solvation free energy of charged molecules based on nonlocal 
electrostatics involving a minimal parameter set which in particular allows to in- 
troduce atomic radii in a consistent way. We test our approach on simple ions and 
small molecules for which both experimental results and other theoretical descrip- 
tions are available for quantitative comparison. We conclude that our approach is 
both physically transparent and quantitatively reliable. 
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1 Introduction 



The unique properties of water, namely its polar nature, high dielectric con- 
stant, and its ability to form hydrogen bonds |lj), are responsible for the 
existence of life as we know it. At the same time, these very properties are 
the main obstacle in modelling hydration effects. An accurate model of ion 
and small molecule hydration will therefore have important applications in 
computational chemistry, chemical engineering, and drug design. 
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In the latter, the interest stems from the fact that drug-hke small molecules 
are usually moderately polar to achieve the required bioavailability. Therefore 
their solvation free energy is generally dominated by the electrostatic contri- 
bution, while the nonpolar contribution is more or less negligible. One of the 
most successful models for predicting solvation free energies of small molecules 
fand proteins) therefore employs the Poisson or Poisson-Boltzmann equations 
Within this approach, the solvent is typically treated as a structureless 
medium of a given dielectric constant e (for water, e ~ 80.) The solute itself 
is treated as a lows "cavity" (0). 

Clearly, this approach fails to take into account many relevant details of 
the solvent structure. The presence of long-ranged electrostatic fields induces 
static correlations among the polar molecules that may vary considerably in 
space depending on the molecular build-up of the solute. These can be taken 
into account by performing Molecular Dynamics simulations, thus leaving the 
continuum description at the expense of high computational cost even for 
small molecules 

Including more structural detail into the continuum description on the one 
hand without increasing computational complexity on the other is therefore 
highly desirable. At the same time, however, also the avoidance of additional 
parameter sets, which is common practice in the use of continuum electrostat- 
ics in computational chemistry, is a key requirement. 

In this work, we follow one particular route to achieve this goal. We develop an 
approach for the computation of the electrostatic (polar) contribution to the 
free energy of solvation based on nonlocal electrostatics (2.; ilO; U; [12) ■ 
This continuum approach, originating in the physics literature as a general- 
ization of classical electrostatics to account for media with spatial dispersion, 
is developed here within a physically transparent minimal parameter set and 
applied to the simplest systems of single ions and small molecules for which ex- 
perimental and theoretical results are available for quantitative comparison. 
In particular we address the issue of a systematic choice of ion radii which 
most commonly is based on essentially empirical parametrizations (0). We 
conclude with an outlook on extensions of our approach for more complex, 
and biologically relevant, molecules. 



2 Nonlocal electrostatics of the solvent 



The electrostatic potential <^ of a charged molecule in solution is, within clas- 
sical (local) electrostatics, given by the Poisson equation (0) 

V [£(r)£oV0(r)] = -e{v) (1) 
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where Eq denotes the dielectric constant of vacuum, and e{r) is a local dielectric 
function which is taken as the constants e^ater ~ 78 of the solution (water), 
and a much smaller value within a solute molecule (for proteins, e.g., a typical 
value taken is e ~ 2 — 5). The rhs of eq.(l) is the charge density giving rise to 
the potential. 

This arguably simple description is prone to complications: first, the transition 
region between solute and solution is ill-described and needs a further (rather 
dehcate) modeling. Second, very little is said about structural effects due to 
the orientation of the polar water molecules near highly charged regions of a 
given solute molecule. 

The simplest way to introduce structural effects into a continuum description 
of the solvent is to account for correlations due to polarization effects between 
solution molecules characterized by a correlation length A. This gives rise to 
a nonlocal generalization of eq.(l) given by (0) 

Vr / dV'e{\r - r'\)eoV,4{r') = -f?(r) (2) 



where isotropy of the liquid medium is supposed. Eq.(2) needs to be accom- 
panied by suitable boundary conditions (see fjlll)). These conditions can be 
simplified considerably by assuming the atomic constituents as conducting 
spheres, which is exact for the case of ions, and approximate for the case of 
small molecules. The solution of the nonlocal solvation problem with a treat- 
ment of the full boundary conditions is beyond the scope of this work and will 



be reported elsewhere fll4 ). 



The rationale behind eq.(2) is the following. Since the electrostatic fields are 
long-ranged, the solvent molecules will feel the presence of the fields of other 
molecules over characteristic distances (that will be of particular relevance 
below). Maintaining locality of the physical fields (i.e., electrostatic potentials 
and electric fields), it is the response functions into which this purely static 
correlation effect needs to be embedded. We note that the dielectric function 
is related to the dielectric response function via e = 1 + x, irrespective of the 
local or nonlocal character of the theoretical description. 

Due to the assumption of spatial isotropy, eq.(2) lends itself to a treatment 
in Fourier space, reducing the linear integro-differential equation eq.(2) to 
an algebraic equation involving the Fourier-transformed function £(|k|) with 
wavevector k. Explicit functional dependences of either e{\r — r'|) or £(|k|) 
have been derived from various approaches, e.g. within a Ginzburg-Landau 
theory for the polarization fields (|lO|). 

Motivated by these works we have formulated a family of nonlocal dielectric 
functions fulfilling the following minimal requirements ([ish : 
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i) for high fc-values (equivalently, on small spatial scales) the solvent molecules 
cannot follow the polarization forces. 



Thus, a limiting value for the macroscopic dielectric function has to be reached 
in this limit. Typical choices in the literature are lim^^oo^l^) ~ ^oo ~ 1-8 
(0; llfih - or = 1 We have considered both cases and do not find a 
significant dependence on this choice; 



ii) in the opposite limit, k —>■ 0, the value for e should equal the macroscopic 
value; 



iii) the physical length scale of the polarization fluctuations is characterized 
by a correlation length A. The local electrostatic limit is given by lim\^oe{\r — 
v'\)=eioc5{r-r'y, 



iv) causality conditions (Kramers-Kronig relations) need to be fulfilled |l7l). 



While the conditions i)-iv) are clearly not sufficient to determine the dielectric 
function unambiguously, we find that the following class of functions in Fourier 
space fulfill these criteria sufficiently well (jl8f ) 



e{k) 



(27r)3/2 



^00 ^lo 



(3) 



which depend on only two parameters, n, and A. For n = 2, the dielectric 
function decays exponentially in real-space, while e{k) is a Lorentzian for 
n = 1, decaying in real-space as exp(— r/A)/r. Note that formally the case 
n = 1 was discussed previously in (j30), based, however, on a different physical 
interpretation (see also (jlfli ) for a similar approach in the context of nuclear 
physics). We have also tested alternative choices, e.g. a Gaussian model and a 
nonlocal model derived from a Ginzburg-Landau theory for polarization fields 
(0). 

The comparison of our nonlocal models with other approaches and experimen- 
tal data is made possible by using the concept of an effective local dielectric 
function e{r). It is usually defined for the potential of a point charge g, 

4)[r) = (4) 
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where 



e{r) 



dk sin(A;r) 



(27r)5/2 kr e{k) 



(5) 



For a general charge distribution e{r) can be defined as well (114)). We note that 
our determination of the effective dielectric constant is consistently performed 
within the nonlocal continuum theory and does not rely on data-based fits 
(SB- 



Expression (5) allows to test our approach against previous results. We took 
empirical models for the radial dependence of e derived for experimental re- 
sults from Mehler and Eichele (ME) (^ and Conway (CO) JH). The CO- 
model is reproduced best by the Lorentzian [n = 1) for A = IsA, while the 
ME-model is reproduced best for A = 24.13A. The purely exponential model 
2) displays unphysical singularities in e{r) at small r, while with a choice 



(n 



of A = 5A, the CO-model is reproduced for larger distances, r > 15A. Figure 1 
shows e{r), as computed by eq.(5) for a single sodium ion, using the Lorentzian 
model with A = IsA. The computation of e{r) performed here makes use of a 
specifically adapted FFT applicable also to molecules ()l4^. 



A Gaussian choice for e{k), however, leads to strong oscillations in the effective 
dielectric function, although its overall shape resembles the CO-model for 
A = SA. Checking our approach against the nonlocal model by Sutmann et al. 
([lol ). we find that their theory fails to comply with our requirement iii), i.e., 
it does not reproduce the correct limiting value at large distances. 



3 Single ions and the choice of their radii 



We now turn to the application of our approach to small ions. We treat these 
ions as Born spheres. A standard problem in the definition of Born-type ions is 
the definition of the ion radius (13). Starting from the accepted interpretation 
that a solvated ion is surrounded by solvation shells, the first of these shells 
will be "as close as possible" to the ion. We can therefore identify the position 
of the centers of the oxygen in the first solvation shell with the first peak in the 
radial distribution function (rdf) which can be obtained either from scattering 
experiments or from molecular dynamics simulations. 

We defined the radius of the water oxygen as half the position of the first peak 
in the oxygen-oxygen radial distribution function (rdf) of bulk water. The ion 
radii were then derived by subtracting this radius from the position of the first 
peak of the ion-oxygen rdf derived from a molecular dynamics simulation. 
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As a specific input for our calculations we take the values for the ion radii 
determined by Aqvist test ). obtained from a combined free eriergy perturba- 
tion/force field approach based on the GROMOS force field (J24'), employing 
several different water models (flexible SPG, rigid SPG, and TIP3P). Based 
on this input, we have determined the solvation free energy from nonlocal 
electrostatics, which is given by 

AG'- = __L/,k{l-^^^^)B^(k) (6) 



where D(k) is the dielectric flux density. 

The result of the calculation is shown in Figs. 2 for monovalent ions and in 
Fig. 3 for divalent ions, for all models to be compared here. Note that we 
chose the parameter A based on its fit to the effective dielectric function. The 
GO-model was reproduced by A = SA for n = 2 (exponential model), by 
A = 15A for n = 1 (Lorentzian) and by A = SA for the Gaussian model. The 
ME-model could only be reproduced by n = 1 with A = 24.13A. The value for 



the exponential model is in accord with the findings in Ref. Also note, 
as shown in Table 1, that the values we obtain are only marginally corrected 
(if at all) by effects due to nonlinear saturation, which contributes in principle 
as well to a reduction of the dielectric constant near the ion. 

The model by Sutmann was used with two different limiting values (Sutmann 
1 with ^oo = 1, and Sutmann 2 with ^oo = 1-8). As the figures show, our 
results are consistently better than all the other theoretical curves. 

The computation of electrostatic solvation free energies is also possible by 
employing eq.(5) in a standard (local) Poisson solver. We have implemented 
this in an available library (BALL) (|25t |lj). Fig. 4 compares the results of 



our computations for the monovalent ions with the nonlocal theory based on 
eq.(6) and from the effective local dielectric function, eq.(5), demonstrating 
the consistency of both approaches. Fig. 5 finally compares the theoretical 
results of the local and nonlocal approaches, obtained with the Poisson solver, 
with the experiment values. Evidently, the nonlocal approach yields results 
consistently superior to those obtained from the standard local theory. 

We stress the significant advantage of our computations to work with first- 
principle radii without arbitrary adjustments. In Figs. 2-5, the Aqvist-radii 
were used (ji^ . In addition, we have also tested the Shannon-radii derived 
from X-ray crystal data (jj?! : lish . without significant effect on our results. 
Our result thus give a basis to the general behef that the need to introduce 
effective radii is in fact a consequence of the local water structure around 



the ions (l26l ). We believe that our nonlocal approach demonstrates that it is 



therefore preferable conceptually to introduce the length-scale governing the 
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structural effects in the solvent rather than introducing necessarily artificial 
procedures to adjust the Born radii. 



4 Small molecules 



As a second application we apply our approach to determine the solvation 
free energy for small alcohols for which we can assume that the polar con- 
tribution exceeds the nonpolar one. Fortunately, accurate measurements of 
the solvation free energies of these molecules are available. We write the 
charge distribution of these molecules as a linear superposition of radially 
symmetric partial distributions translated by a vector R, Qi = f)(|r + Rj|), i.e. 
f?(r) = J2iLi ft(|r + RiD- Its Fourier transform is then given by the expression 
^'(k) = Eili Qi{k) exp(-iRi ■ k). 

Again we have to define the radii of the atoms. For this we chose to classify the 
atoms into classes depending on their chemical environment in the molecule 
(e.g., like the hydrogen atom in an OH-group), expecting that the radii of all 
atoms in a certain class are more or less similar. For methanol (CH^OH), e.g., 
we defined four classes of atoms and used their rdf 's with the oxygen of water 
to find the following set of radii: Methyl C: 2.135 A, hydroxyl O: 2.014 A, 
hydroxyl H: 1.115 A, methyl H: 1.394 A. 

Compared to the calculation of the free energies for the ions this computation 
is slightly more involved as it requires a three-dimensional integration instead 
of a one-dimensional one. Within the framework of the simplified boundary 
conditions, the unperturbed fields of the atomic constituents are simply su- 
perimposed. To ensure short computation times, in our implementation we 
have used the VEGAS-Monte-Carlo i nteg ration scheme (0) that is supplied 
with the GNU scientific library GSL (|30l). Details of these computations will 



again be given elsewhere ()l4^. 



The results for some small alcohols are shown in Table 1. For the comparison 
we have used the exponential with A = 5 A, the Lorentzian with A = 15 A, 
and the Lorentzian with A = 24.13 A. All values in these tables are given in 
kJ/mol. 

The interpretation of the results is more complicated than in the Born ion 
case. While there we could assume that the contribution of the nonpolar part 
of the free energy of solvation could be neglected (for the Born ions it can be 
estimated to be of the order of 10 — 20 kJ/mol), this is not the case for the 
alcohols. The electrostatic contribution is still the dominant part, but not the 
only significant one. We therefore applied a very simple model for the nonpolar 



contribution (|3l[ ) in order to be able to compare our results to experimental 
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data. The results given for the polar contribution can also be optimized by 
considering improved charge distributions (14). 



5 Conclusions and outlook 

We have demonstrated that the approach to the computation of solvation free 
energies based on nonlocal electrostatics is able to reproduce experimental 
data for ions and small molecules with reasonable accuracy. The approach we 
put forward has the basic advantage to rely on essentially one parameter which 
has a transparent physical interpretation as the correlation length of the po- 
larization fluctuations. While within the present paper this parameter is used 
as a fitting parameter, it should be clear that it might be also determined by 
experiment or simulation. A further significant advantage is that our approach 
does not make use of the commonly used adjustments of atomic or ionic radii. 
A challenge for the future will be the extension of our approach to determine 
electrostatic solvation free energies for more complex and biologically relevant 
molecules. Work in this direction is in progress (|l4l). 

Acknowledgement. This work is supported by the DFG under Schwerpunkt- 
sprogramm "Informatikmethoden zur Analyse und Interpretation grosser geno- 
mischer Datenmengen" (grant LE952/2-1). 



References 

[I] P. Ball, H2O: A biography of water, Weidenfeld and Nicholson (2000) 
[2] B. Honig, K. Sharp, A.-S. Yang, J. Phys. Chem. 97, 1101 (1993) 

[3] B. Honig, A. Nichohs, Science 268, 1144 (1995) 
[4] R.M. Jackson, M.J.E. Sternberg, J. Mol. Biol. 250, 258 (1995) 
[5] R.M. Levy, E. Gallicchio, Annu. Rev. Chem. 49, 531 (1998) 
[6] J.D. Jackson, Classical Electrodynamics, John Wiley & Sons, Inc. 3rd ed. 
(1998) 

[7] R.R. Dogonadze, A. A. Kornyshev, A.M. Kuznetsov, Teor. Mat. Fiz. 15, 
127 (1973) 

[8] A.A. Kornyshev, A.G. Volkov, J. Electroanal. Chem. 180, 363 (1984) 
[9] A.A. Kornyshev, in The Chemical Physics of Solvation, R.G. Dagodnaze 

et al. (eds.), Elsevier, Amsterdam (1985) 
[10] G. Sutmann, Die nichtlokale dielektrische Funktion von Wasser, PhD - 

Thesis, Jiilich (1999) 

[II] M.V. Basilevsky, D.F. Parsons, J. Chem. Phys. 108, 9107 (1998) 
[12] M.V. Basilevsky, D.F. Parsons, J. Chem. Phys. 108, 9114 (1998) 



8 



[13] M.S. Lee, F.R. Salsbury Jr., C.L. Brooks III, J. Chem. Phys. 116, 10606 
(2002) 

[14] A. Hildebrandt, O. Kohlbacher, R. Blossey, H.P. Lenhof, in preparation 
(2003) 

[15] A. Hildebrandt, An algorithm for the efficient and reliable computation of 
the electrostatic contribution to the free energy of solvation using nonlocal 
electrodynamics, Diploma thesis, Saarland University (2002) (unpublished) 
[16] J.B. Hasted, Aequous Dielectrics, Chapman & Hall, London (1973) 
[17] L.D. Landau, E.M. Lifschitz, Electrodynamics of Continuous Media, Perg- 

amon Oxford (1960) 
[18] We have normalized the Fourier transform symmetrically. 
[19] U. Ritschel, L. Wilets, J.J. Rehr, M. Grabiak, J. Phys. G: Nucl. Part. 

Phys. 18, 1889 (1992) 
[20] E.L. Mehler, G. Eichele, Biochemistry 23, 3887 (1984) 
[21] B. Mallik, A. Masunov, T. Lazaridis, J. Comp. Chem. 23, 1090 (2002) 
[22] B.E. Conway, J.O.M. Bockris, LA. Ammar, Trans. Faraday Soc. 47, 756 
(1951) 

[23] J. Aqvist, J.Phys. Chem. 94, 8021 (1990) 

[24] W.F. van Gunsteren and H.J.C. Berendsen, Groningen Molecular Simu- 
lation (GROMOS) Library Manual Version 3.0 

[25] O. Kohlbacher, H.-P. Lenhof, Bioinformatics 16, 815 (2000) 

[26] L. Sandberg, O. Edholm, J. Chem. Phys. 116, 2936 (2002) 

[27] R.D. Shannon, C.T. Prewitt, Acta Crystallogr., Sect. B: Struct. Crystal- 
logr. Cryst. Chem. 25, 925 (1969) 

[28] R.D. Shannon, Acta Crystallogr. Sect. A: Cryst. Phys., Diffr., Theor. 
Gen. Crsytallogr. 32, 751 (1976) 

[29] G.P. Lepage, J. Comp. Phys. 25, 192 (1 978) 

[30] to be found at: http:/ /sources. redhat.com/gsl| 

[31] H. Uhlig, J. Phys. Chem. 41, 1215 (1937) 

[32] Y. Marcus, Ion Properties, Dekker, New York (1997) 



9 




5 10 15 20 25 

r [Angstrom] 



Fig. 1. e(r) for a single sodium ion 
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Results from Eq. (6) [kJ/mol] 

Fig. 4. Free energy of solvation for monovalent ions, from the nonlocal expression, 
eq.(6), and from the solution of the local Poisson equation with the effective dielec- 
tric constant, eq.(5). 




Fig. 5. Free energy of solvation for monovalent ions, from the nonlocal expression, 
eq.(6), (circles) and from the solution of the local Poisson equation (squares), i.e. 
the Born energies. (The experimental values corresponding to the ions are on the 
diagonal.) 
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Ion 


r[A] 


AGs 






corr 




calc 


^Gexp 


Li+ 


0.645 


-1063.08 


-568.9 


5.03 





-563.87 


-563.87 


-481 


Na+ 


1.005 


-682.30 


-392.51 


5.48 





-387.03 


-387.03 


-375 


K+ 


1.365 


-502.36 


-306.79 


6.01 





-300.78 


-300.78 


-304 


Rb+ 


1.505 


-455.63 


-284.02 


6.23 





-277.79 


-277.79 


-281 


Cs+ 


1.715 


-399.84 


-256.39 


6.59 





-249.8 


-249.8 


-258 


Mg2+ 


0.615 


-4459.74 


-2370.66 


4.99 


110 


-2365.67 


-2255.67 


-1838 


Ca?+ 


1.015 


-2707.32 


-1557.41 


5.49 


45 


-1551.92 


-1506.92 


-1515 


Sr2+ 


1.195 


-2295.29 


-1364.6 


5.75 


30 


-1358.85 


-1328.85 


-1386 


Ba2+ 


1.385 


-1980.43 


-1213.12 


6.04 


24 


-1207.08 


-1183.08 


-1259 



Table 1 

Comparison of different models for the hydration free energy of different mono- 
and divalent cations. All energies are in kJ/mol. AGb is the hydration free energy 
computed for a Born ion of radius r. AGnioc is the electrostatic hydration free 
energy computed with our nonlocal Lorentzian model (A = 24.13). AGnpoi is the 
nonpolar contribution to the hydration free energy as computed with the model 
of Uhlig (ilj). AG 

corr is the nonlin ear correction to the hydration free energy for 
the corresponding ion (taken from (|2^). AGexp is the experimental hydration free 
energy (taken from (Q). AGcaic = AG nioc + AGnpoi- AG^°[^^ = AGcaic + ^^Gcorr- 
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Substance 


Exponential 
A = 5 A 


Nonpolar 


E 


Experiment 


Ethanol 


-32.84 


8.89 


-23.95 


-20.51 


Methanol 


-25.64 


8.10 


-17.54 


-21.26 


1-Butanol 


-29.79 


10.31 


-19.48 


-19.76 


1-Hexanol 


-37.41 


11.83 


-25.58 


-18.25 


Oct auol 


-33. (33 


13.28 


-20.35 


-17.12 


Cyclopentanol 


-28.90 


10.31 


-18.59 


-22.98 




Lorentzian 
A = 15 A 








Ethanol 


-34.17 


8.89 


-25.28 


-20.51 


Methanol 


-29.49 


8.10 


-21.39 


-21.26 


l-Bu(anol 


-31 13 


10.31 


-20.82 


-19.7G 


1-HexanoI 


-38.88 


11.83 


-27.05 


-18.25 


Octanol 


-35.08 


13.28 


-21.79 


-17.12 


Cyclopentanol 


-30.16 


10.31 


-19.85 


-22.98 




Lorentzian 
A = 24.13 A 








ElllMllol 


-32.89 


8.89 


-21.0 


-20.51 


Methanol 


-27.88 


8.10 


-19.78 


-21.26 


1-Butanol 


-29.94 


10.31 


-19.63 


-19.76 


1-Hexanol 


-37.65 


11.83 


-25.82 


-18.25 


Octanol 


-33.77 


13.28 


-20.49 


-17.12 


Cyclopentanol 


-29.06 


10.31 


-18.75 


-22.98 



Table 2 

Results for the free energy of solvation for small alcohols. All free energies are given 
in kJ/mol. 
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